library(RSocrata)
library(tidyverse)
library(lubridate)
library(tidygeocoder)
library(rgeoboundaries)
library(sf)
library(tmap)
library(spatstat)
source("Funciones_propias_ppp.R", encoding = "UTF-8")
La accidentalidad es uno de los problemas más comunes en el planeta provocando perdidas que pueden ir desde materiales hasta la vida misma, por lo que es de particular interés tratar de comprender como se comporta este fenómeno para evitar ser victimas de él.
La problemática en cuestión tiene la complicación de ser algo de naturaleza aleatoria, sin embargo esto no quiere decir que no exista alguna forma de controlarla pues la experiencia ha mostrado que las diferentes ciudades tienen lugares con mayor concentración de accidentes que otras ubicaciones dentro de la misma ciudad.
Debido a que se está trabajando la ocurrencia de accidentes en una ciudad de interés, la distribución Poisson y los procesos puntuales Poisson resultan ser una herramienta adecuada para tratar de explicar el fenómeno.
El estudio de accidentalidad se hará en los años 2019, 2020 y 2021 con el propósito de ver como se comporta la accidentalidad prepandemia, pandemia y postpandemia.
Una vez presentado el fenómeno que se trabajará, se selecciona la ciudad de Barranquilla para realizar el estudio.
Los datos se extrajeron de la página web GOV.CO la cual tiene diferentes bases da datos de Colombia abiertas al público. ( Accidentalidad en Barranquilla). La base de datos contiene una gran cantidad de variables, sin embargo solo se escogen aquellas de interés para el patrón puntual. La estructura de los datos es presentada a continuación:
accidentes <- read.socrata("https://www.datos.gov.co/resource/yb9r-2dsi.json") %>%
mutate(fecha = ymd(fecha_accidente)) %>%
select(12, 6:8) %>%
filter(year(fecha) > 2018 & year(fecha) < 2022)
knitr::kable(head(accidentes, 10),
col.names = c("Fecha", "Gravedad del accidentes",
"Clase de accidentes", "Sitio accidente"))
| Fecha | Gravedad del accidentes | Clase de accidentes | Sitio accidente |
|---|---|---|---|
| 2019-01-01 | Solo daños | Choque | CALLE 82 CRA 71 |
| 2019-01-01 | Con heridos | Choque | CL 30 CR 27 |
| 2019-01-01 | Con heridos | Choque | CL 17 CR 30 |
| 2019-01-01 | Solo daños | Choque | CRA 31 CALLE 68C |
| 2019-01-01 | Con heridos | Atropello | CL 99D CR 8E |
| 2019-01-01 | Solo daños | Choque | VIA 40 CR 67B |
| 2019-01-02 | Solo daños | Choque | AV CIRCUNVALAR FRENTE AL NO 70A 12 |
| 2019-01-02 | Con muertos | Choque | AV CIRCUNVALAR CR 13 |
| 2019-01-02 | Solo daños | Choque | AV CIRCUNVALAR 49 39 |
| 2019-01-02 | Solo daños | Choque | CR 53 68 216 |
Se hace necesario realizarle unos ajustes a la base de datos
considerada puesto que está no incluye las ubicaciones exactas de los
accidentes ya sea en longitud - latitud o proyección UTM, para dicho
propósito se usa la función geo() del paquete
tidygeocoder para convertir las direcciones de los
accidentes en coordenadas de longitud - latitud.
direcciones <- paste(accidentes$sitio_exacto_accidente, ", Barranquilla, Colombia", sep = "")
localizaciones <- geo(address = direcciones, method = "arcgis")
write.csv(x = localizaciones, file = "accidentes.csv", row.names = F)
Luego de realizar obtener las coordenadas en longitud - latitud, se
usan múltiples funciones del paquete sf para obtener las
respectivas ubicaciones en formato UTM obteniendose lo siguiente:
bqlla <- geoboundaries(country = "COLOMBIA", adm_lvl = 2) %>%
filter(shapeName == "DISTRITO ESPECIAL, INDUSTRIAL Y PORTUARIO DE BARR*")%>%
st_transform(crs = 3857)
coordenadas <- read.csv("accidentes.csv")
coordenadas <- cbind(accidentes, coordenadas[, 2:3]) %>%
st_as_sf(coords = c('long', 'lat')) %>%
st_set_crs(value = 4326) %>%
st_transform(crs = 3857) %>%
st_intersection(bqlla)
accidentes <- list()
for(i in 2019:2021){
accidentes[[as.character(i)]] <- filter(coordenadas, year(fecha) == i & month(fecha) > 3)
}
knitr::kable(head(accidentes[["2019"]][, c(1:4, 10)], 10))
| fecha | gravedad_accidente | clase_accidente | sitio_exacto_accidente | geometry |
|---|---|---|---|---|
| 2019-04-01 | Solo daños | Choque | CL 56 14 30 | POINT (-8327234 1227521) |
| 2019-04-01 | Con heridos | Choque | CL 30 CR 25 | POINT (-8324772 1228159) |
| 2019-04-01 | Con heridos | Choque | CR 27 CL 85 | POINT (-8330025 1229473) |
| 2019-04-01 | Solo daños | Choque | CR 38 CL 81 | POINT (-8329330 1230390) |
| 2019-04-01 | Solo daños | Choque | CALLE 31 CRA 16 | POINT (-8324126 1229674) |
| 2019-04-01 | Solo daños | Choque | CL 90 90B 42 | POINT (-8329448 1226070) |
| 2019-04-01 | Solo daños | Choque | CL 94 CR 71 | POINT (-8328610 1234554) |
| 2019-04-01 | Con heridos | Choque | AV CIRCUNVALAR CL 74 | POINT (-8331175 1231848) |
| 2019-04-01 | Con heridos | Choque | AV CIRCUNVALAR CR 13 | POINT (-8331175 1231848) |
| 2019-04-01 | Solo daños | Choque | CL 40 CR 24 | POINT (-8325400 1228241) |
tmap_mode('view')
tm_shape(bqlla)+
tm_polygons(alpha = 0.3, border.alpha = 0.7)+
tm_shape(shp = accidentes[['2019']])+
tm_dots(size = 0.01)
tmap_mode('view')
tm_shape(bqlla)+
tm_polygons(alpha = 0.3, border.alpha = 0.7)+
tm_shape(shp = accidentes[['2020']])+
tm_dots(size = 0.01)
tmap_mode('view')
tm_shape(bqlla)+
tm_polygons(alpha = 0.3, border.alpha = 0.7)+
tm_shape(shp = accidentes[['2021']])+
tm_dots(size = 0.01)
Uno de los parámetros críticos al momento de modelar un patrón puntual es la intensidad la cual puede ser constante (homogénea) o variable (inhomogénea), por lo tanto es importante iniciar el análisis con una prueba de homogeneidad de la intensidad.
# Guardando datos por anio
datos_2019 <- ppp(x = st_coordinates(accidentes[[1]])[, 1],
y = st_coordinates(accidentes[[1]])[, 2],
window = as.owin(W = bqlla))
datos_2020 <- ppp(x = st_coordinates(accidentes[[2]])[, 1],
y = st_coordinates(accidentes[[2]])[, 2],
window = as.owin(W = bqlla))
datos_2021 <- ppp(x = st_coordinates(accidentes[[3]])[, 1],
y = st_coordinates(accidentes[[3]])[, 2],
window = as.owin(W = bqlla))
Para iniciar, se divide el mapa de Barranquilla en 25 cuadrantes, si el patrón fuera homogéneo se esperaría que el número de accidentes en cada uno de las subsecciones de la ciudad contenga aproximadamente la misma cantidad de accidentes.
qc_2019 <- quadratcount(datos_2019, nx = 5, ny = 5)
plot(qc_2019, main = "Accidentes en 2019")
qc_2020 <- quadratcount(datos_2020, nx = 5, ny = 5)
plot(qc_2020, main = "Accidentes en 2020")
qc_2021 <- quadratcount(datos_2021, nx = 5, ny = 5)
plot(qc_2021, main = "Accidentes en 2021")
En todos los años se aprecia que hay una gran discrepancia entre el número de accidentes por sectores. Claramente al este de la ciudad el número de accidentes es mayor respecto al oeste lo cual es una fuerte señal de que el patrón puntual en cuestión es de naturaleza inhomogénea.
Adicional a los conteos del número de accidentes en cada uno de los sectores definidos previamente, se usa la prubea \(\chi^2\) para contrastar:
\[ \begin{cases} H_0: \text{La intensidad es constante} \\ H_1: \text{La intensidad no es constante} \end{cases} \] A continuación se presentan los resultados obtenidos
quadrat.test(qc_2019)
##
## Chi-squared test of CSR using quadrat counts
##
## data:
## X2 = 2857.1, df = 18, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 19 tiles (irregular windows)
quadrat.test(qc_2020)
##
## Chi-squared test of CSR using quadrat counts
##
## data:
## X2 = 1324.7, df = 18, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 19 tiles (irregular windows)
quadrat.test(qc_2021)
##
## Chi-squared test of CSR using quadrat counts
##
## data:
## X2 = 2094.2, df = 18, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 19 tiles (irregular windows)
Puesto que en todos los casos se tienen valores - p demasiado pequeños (orden de \(10^{-16}\)) se rechaza la hipótesis de homogeneidad de intensidad en todos los casos, por lo tanto se hace necesario estimarla con algún método.
A continuación se presentan mapas de probabilidad de accidentalidad en los años considerados al realizar una estimación no parametrica de la función de intensidad basado en funciones kernel.
La selección del ancho de banda se puede escoger de muchas maneras y existen multiples métodos para calcularlos. Luego de usar diferentes propuestas para el calculo de este y probar con valores empíricos se llego a que un valor razonable para el mismo es 350 pues captura la concentración de accidentes de forma plausible.
graph_ppp(datos_2019, 350, "Año 2019")
graph_ppp(datos_2020, 350, "Año 2020")
graph_ppp(datos_2021, 350, "Año 2021")
Previamente se verifico que el número de accidentes en la ciudad tiene intensidad inhomogénea, en esta sección el propósito es estimar dicha función usando modelos log-lineales los cuales son de la forma
\[ log \ \lambda_\theta (u) = \theta \cdot S(u) \]
donde \(S(u)\) es una función vectorial de valor real evaluada en alguna ubicación \(u\).
Para ajustar estos modelos en R, se usan la función
ppm() del paquete spatstat la cual recibe como
primer argumento un objeto de la clase ppp y como segundo
una formula.
El ajuste de dichos es modelos es sensible a la escala de medición de las ubicaciones; debido a que las ubicaciones se encuentran en metros (los cuales son números muy grandes) es necesario cambiar la escala de medición por una más “pequeña” para evitar problemas de singularidad matricial por lo que se convierten los metros en kilometros.
# Rescalando los objetos ppp para ajustar los modelos
# Ano 2019
datos_2019 <- rescale(datos_2019, 1000)
unitname(datos_2019) <- c("km", "kilometers")
# Ano 2020
datos_2020 <- rescale(datos_2020, 1000)
unitname(datos_2020) <- c("km", "kilometers")
# Ano 2021
datos_2021 <- rescale(datos_2021, 1000)
unitname(datos_2021) <- c("km", "kilometers")
Luego de reescalar las ubicaciones se ajustan algunos modelos log-lineales.
# Ajuste de modelos
# Tendencia lineal en x e y
log_lin_2019 <- ppm(datos_2019, ~ x + y)
log_lin_2020 <- ppm(datos_2020, ~ x + y)
log_lin_2021 <- ppm(datos_2021, ~ x + y)